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I. INTRODUCTION 



Spiral waves are a type of self-organization observed in 
a large variety of spatially extended, thcrmodynamically 
non-equilibrium systems of physical, chemical and biolog- 
ical nature [1-19], where wave propagation is supported 
by a source of energy stored in the medium. If the system 
can be considered spatially uniform and isotropic and 
its properties do not depend on time, the corresponding 
mathematical models possess corresponding symmetries. 
For many practical applications, a considerable interest is 
in non-stationary dynamics of spiral waves, which is usu- 
ally defined separately either as drift, which is displace- 
ment of the average position of the core of the spiral with 
time due to external symmetry-breaking perturbations, 
or meandering, which is spontaneous symmetry break- 
ing due to internal instability rather than external forces 
and which is manifested by complicated movement of the 
spiral with the average position of the core typically un- 
moved. 

The numerical simulation of drift and meander of spi- 
ral waves, particularly when models are complicated and 
high accuracy is required, can be challenging. There are 
some theoretical considerations which suggest some way 
of dealing with this challenge. So it has been observed 
that as far as drift is concerned, spiral waves behave like 
particle-like objects, which results from effective localiza- 
tion of the critical cigenfunctions of the adjoint linearized 
operator [20-24] , so it should be sufficient to do the com- 
putations only around the core of the spiral to predict 
its drift. On the other hand, in the absence of external 
symmetry breaking perturbations, meandering of spirals 
can be understood by explicitly referring to the Euclidean 
symmetry of the unperturbed problem [25-30] . Specifi- 
cally, an idea of dynamics in the space of symmetry group 
orbits [31], when applied to a reaction-diffusion system 
of equations and the Euclidean symmetry group, leads to 
a description which is formally equivalent to considering 
the solution in a moving frame of reference (FoR) such 
that the spiral wave maintains a certain position and ori- 
entation in this frame [29] . We shall call it comoving FoR 



for short. 

The purpose of this article is to present a computa- 
tional approach based on these considerations. We cal- 
culate the dynamics of the spiral wave in a comoving FoR; 
as a result, the core of the spiral never approaches the 
boundaries of the computation box, which allows compu- 
tations of drift and meandering of large spatial extent us- 
ing small numerical grids. A simple software implementa- 
tion of this approach, which is based on the popular spi- 
ral wave simulator 'EZ-SPIRAL' [32, 33], and which we 
called 'EZRide', is provided on the authors' website [34]. 

Our approach can be compared to the approach pro- 
posed by Beyn and Thummlcr [35] and further developed 
by Hermann and Gottwald [36]. Their approach also ex- 
ploits symmetry group orbits, but is different in some 
essential details. We shall discuss the similarities and 
differences when we will have introduced our method. 

The structure of the paper is as follows. In section II we 
lay out mathematical basics of the approach and briefly 
compare it with [35]. In section III we describe the nu- 
merical method itself. In section IV we illustrate the work 
of the method by simple and quick examples. The poten- 
tial for numerical accuracy is demonstrated in section V. 
The subsequent three section are dedicated to examples 
of applications of the methods to problems where the 
conventional methods would be struggling: section VI 
for the degenerate case of meandering which results in 
"spontaneous drift" of spirals; section VII for the dy- 
namics near to, and beyond, the parametric boundary 
at which the core radius of the spiral becomes infinite; 
and section VIII for drift caused by a symmetry break- 
ing perturbation. We conclude with a brief discussion of 
the results in section IX. 



II. SYMMETRY GROUP REDUCTION 

Following [29], we start from a perturbed reaction- 
diffusion system of equations in a plane, 



du 

~dt 



DV 2 u + f (u) + h(u, Vu, r, t), 



(1) 
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where u = (u^\ . . . , u^) = u(r, t) G R™ is a column- 
vector of reagent concentrations varying in space and 
time, f = f(u) is a column- vector of reaction rates, 
D G R nx ™ is the matrix of diffusion coefficients, h G R™ 
represents symmetry-breaking perturbations, ||h|| <C 1, 
n > 2, and r = (a;, y) G K 2 . If h = 0, then equation (1) 
is equivariant with respect to Euclidean transformations 
of the spatial coordinates r. 

The following technical discussion is necessary to place 
our method in the context of other works in the field. 
Readers not interested in technical details may skip down 
to system (16). 

The idea of the symmetry group reduction is conve- 
nient to describe if we view (1) as an ordinary differential 
equation in a suitably chosen functional space B, 

^H=F(U) + H(U,t) (2) 

where U : R — > B represents the dynamic field u, F : B — > 
B represents the unperturbed right-hand side DV 2 u + f , 
and H : B x R — >• B represents the perturbation h. 

Let us suppose that equation (2) at h = is equivari- 
ant with respect to a representation T of a Lie group Q 
in B. This means that for any g G Q and any U € B, we 
have 

F(T( 5 )U) = T( 5 )F(U). (3) 

In our case, Q = SE(2), the special Euclidean transfor- 
mations of the plane R 2 — > R 2 (including translations 
and rotations), and T is its representation in the space 
of functions u(r) defined on this plane, acting as 

T(.g)u(r) = uGTV). (4) 

We consider a subset Bq C B such that Q acts freely on 
Bo, i.e. for a U € Bq, any nontrivial transformation g G Q 
changes U, in other words, T(<?)U = U =>■ g = id. In 
the terminology of [31], Bq is the principal stratum of B, 
corresponding to the trivial isotropy subgroups. In our 
case, this means that the graph of the function u(r) £ 
Bq is devoid of any rotational or translational symmetry, 
which is of course true for functions describing single- 
armed spiral waves. 

It is straightforward that at H = 0, the set Bq is an 
invariant set of (2). Moreover, we shall restrict our con- 
sideration to such perturbations H(f) that resulting so- 
lutions U(i) remain in Bo for all t. This means, that the 
perturbations are supposed to be so small they cannot 
impose incidental symmetry on the otherwise unsymmet- 
ric spiral wave solutions. 

A group orbit of a given U is defined as the set 
T(£)U = {T(g)V | g G Q}. That is, it is a set of all such 
functions u(r) that can be obtained from one another 
by applying an appropriate Euclidean transformation to 
r. A group orbit is a manifold in Bq, of a dimensionality 
equal to d — dim Q less the dimensionality of the isotropy 
group. In our case, dimSE(2) = 3, the isotropy group is 



trivial and the orbits are smooth three-dimensional man- 
ifolds. 

From the definition of the set Bq it follows that this 
set is foliated by group orbits. The principal assumption 
for the following analysis is that there exists an open 
subset S C Bq, also invariant with respect to Q, in which 
the foliation has a global transversal section, i.e. we can 
select one representative from each orbit in S, such that 
all such representatives form a smooth manifold A4 C S, 
which is everywhere transversal to the group orbits. We 
call this manifold a Representative Manifold (RM). That 
would mean that any orbit in S crosses M. transversally 
and exactly once. Hence 

VUgS, 3'(g,V) G G x M : U = T(<?)V. (5) 

The RM has co-dimensionality equal to the dimensional- 
ity of the group orbits, i.e. in our case codim A4 = d = 3. 
It is assumed to be smooth and we expect that it can 
locally be described by equations /Xf (V) = 0, £ = 1, . . . d, 
where functions fi# : B — > R, i. e. are functionals when in- 
terpreted in terms of the original reaction-diffusion equa- 
tion (1). 

A convenient pictorial interpretation for our case is in 
terms of spiral wave solutions and their tips. Suppose the 
conditions /Ui(V) = 0, /U2(V) = determine that the tip 
of the spiral wave is located at the origin, and condition 
/i 3 (V) = fixes its orientation, so M consists of such 
functions that look like spiral waves V which have the tip 
exactly at the origin and in a standard orientation. Then 
equation (5) states that any spiral wave solution u(r), 
considered at a fixed moment of time, can be transformed 
by a Euclidean transformation, in a unique way, to a 
solution v(f) which has its tip at the origin and in the 
standard orientation. This is equivalent to saying that 
v(r) is the same as u(r) only considered in a different 
system of coordinates, with the origin at the tip of u(r) 
and oriented accordingly to the orientation of that tip. 
We shall say this is the system of coordinates attached to 
the tip. An example of \xi, as used e.g. in [29], is 

A*i[v(0]=« (,l) (3)-'u., (H 
A*a[v(» ! )]=» (,a) (S)-»„ (6b) 
M 3 [v(f)] =d x v^(0), (6c) 

where {li,h,h} C {l,...,n} are suitably chosen com- 
ponents, and l\ 7^ 12- This means that the tip of u(r) 
is defined as the point of intersection of isolines of the 
components l\ and I2 of the field u at appropriately cho- 
sen levels it* and u* respectively, and the orientation of 
the attached coordinate system is such that gradient of 
component I3 (which may or may not coincide with l\ or 
I2) is along the y-axis in that system. This choice of lit 
is of course not prescriptive, and later in this paper we 
shall consider some variations. 

Regardless of the exact definition of the tip, i.e. choice 
of functionals [in, an essential assumption that we have to 
make is that our spiral waves have one tip only, otherwise 
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sides, we get 



FIG. 1: (Color online) Sketch of skew-product decomposi- 
tion of an equivariant flow using a Representative Manifold 
M, which has exactly one transversal intersection with evey 
group orbit g £ Q within the relevant stratum of the phase 
space B and is diffeomorphic to the orbit manifold. Trajectory 
(U,U',U") of an equivariant flow in B is a relative periodic 
orbit, since it projects onto the trajectory (V,V',V" = V) 
on M which is periodic. The flow on M is devoid of symmetry 

g. 



there would be more than one way to transform them to 
the standard position or to chose the attached system of 
coordinates. Hence the reason for a further constraint to 
the subset S C Bo . which we now can define as consisting 
of such one-tip spiral wave solutions, or functions that 
look like it: without such constraint, the whole set Bo 
includes solutions with no tips or more than one tip, for 
which the decomposition (5) would not hold. As before, 
we assume that set S is invariant with respect to the 
dynamic equation (2) for not too big ||H||, that is, if 
U(0) G S, then U(<) g S for all t > and ||H|| < # max . 

A further restriction is on the manifold M. It is easy 
to see that equations like (6) may not be sufficient to 
define the manifold with the required property that any 
orbit crosses it only once. For instance, if v(r) satisfies 
(6), then v(— r) also satisfies it, so a rotation by 180° 
around the origin transfers a point on S to another point 
on S. So to make the representation (5) unique, rather 
than just requiring that the gradient of the ^-component 
of v(r) is along the y axis, one would need to specify in 
which direction it is, say add to the definition of M. by 
the equations ^1,2,3 [v] = a further inequality 

/i 4 [v]>0, where /i 4 [v(r)] = d y v ih) . (6d) 

This comment extends to the variations of (6) which we 
consider later. 

By performing the decomposition (5) for every t > 0, 
we decompose motion in S to motion along the RM and 
motion along group orbits which are transversal to the 
RM (see illustration in fig. 1). 

So for alH > 0, we have 



U(t) = T( ff )V(t) 



(7) 



T( ^ 1} v + tSf = F(v) + fi(v ' 9 ' t] (8) 



where 



H(V, ff ,f) =T(g- 1 )U(T(g)V,t) 



(9) 



We note that if H = 0, the right-hand side of (8) is 
independent of g. 

By the assumptions made, intersection of the group 
orbit T(C?)V with the manifold M. at the point V is 
transversal. This means that the vectors F(V) and 
H(V, g, t) can be uniquely decomposed into the sums of 
the components along the group and along the manifold, 

F(V) = F g (V)+F M (V), (10a) 
H(V, 5 ,t)=H (V > t)+H M (V > t). (10b) 



Hence equation (8) splits into two components, along 
the RM and along the group orbit (GO): 



(RM) 



<9V 
~dt 



F M (v) + n M (v,t) 



(lla) 



Substituting (7) into (2) and applying T (g x ) to both 



(GO) T( 5 - 1 )^Mv = F e (V)+H g (V,t) (lib) 

Note that equation (lla) is the equation on the infinite- 
dimensional manifold yVf, i.e. corresponds to a partial 
differential equation, whereas the left- and right-hand 
sides of the equation (lib) are in the tangent space to the 
finite-dimensional group orbits, and the dynamic variable 
g is an element of the finite-dimensional manifold Q, so 
(lib) is in fact a system of ordinary differential equations 
of order d = dim Q . 

At this point we comment on what we see as a signifi- 
cant difference between our approach and that proposed 
by Beyn and Thummler [35] (BT). Using our notation, in 
place of our "pinning" conditions ne(V) = 0, I = 1, . . . d, 
they defined "phase conditions" of the form £tf(V, 5) = 
(see equation (2.22) in [35]), subsequently further gen- 
eralized to fie(V,g,dg/dt) = (ibid., equation (2.33)). 
This means that their decomposition U = T(g)V is not 
uniquely determined by the current state U, but depends 
on history. Such generalization may have its advantages 
and, apparently, works well for relative equilibria, i.e. 
steadily rotating spirals [35, 36]. However, the situation 
is different if the solution is a meandering spiral, i.e. is 
periodic with period P in the orbit space (as illustrated 
in fig. 1). This means that U(i + P) is equivalent to 
U(t) up to some Euclidean transformation. In our ap- 
proach, it is then guaranteed, that V(< + P) = V(f), 
as by (5), fie(T(g^ 1 )\J) = has a unique solution for g 
at a given U. However, in the BT approach, typically 
V(t + P) ^ V(t), since ^(Tfo-^U, g, dg/dt) = docs 
not uniquely define g, as dg/dt is not fixed. So in our 
approach, study of meandering spirals reduces to study 
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of periodic solutions for V(i), but it does not do so in 
the BT approach. 

A practical approach to the problem of decomposing 
the vector fields as in (10) is as follows. Equations (11) 
together with the definition of the RM via functionals \xi 
can be re-written in an equivalent form 



H 



dV _ 
~dt ~ 

w(v(t))=o, e = i,.. 



,d, 



(12a) 
(12b) 

(12c) 



where A = A(V,t) = -F e (V) - H e (V,i) is a vector 
belonging to the three-dimensional tangent space of the 
group orbit T(Q)(V) at V. In this formulation, at any 
given moment of time, equations (12a) and (12b) to- 
gether define the evolution of V and the current value 
of the vector A, whereas equation (12c) defines the evo- 
lution of g. 

By definition, vector A is a result of action of a linear 
combination of the generators of the Lie group T(G) as 
linear operators on V. To write the explicit expression 
for the general form A for our case, let us introduce co- 
ordinates (R, 9) on g = SE{2), where R = {X, Y) is the 
translation vector, 6 is the rotation angle and a group 
clement acts as 



g = (R,Q) : r^R 



(13) 



where 7 



-1 

1 



, so exp(70) is the matrix of rotation 



by angle 0. 

Using this representation, differentiating the definition 
of T(g)v given by (4), and substituting the result into 
(12c), we get 



where 



cudev + (c ■ V)v, 



= 6, c = e^ e R 



(14) 



(15) 



and 8 is the polar angle in the (x, y) plane, so de = 
xd y - yd x . 

With this result, the system (12) in the original PDE 
notation states 



~dt 



= DV 2 v + f (v) + h ( v, e^Vv, R + e 7t V, * 



dx 



de 



0. 

dR 
dt ~~ 



V)v + c-, 


(16a) 


w ( ' 2) (0,t) = u», 


(16b) 


dy 


(16c) 




(16d) 



where the dynamic variables are v(r, i), c(t), u>(t), R{t) 
and Q(t). 

In terms of the tip of the wave, equation (16a) is the 
original reaction-diffusion equation (1) written in the co- 
moving FoR, equations (16b), (16c) define the attach- 
ment (pinning) of the tip to this FoR, and equations (16d) 
describe the movement of the FoR and, therefore, of the 
tip. 

Equations (16b), (16c) imply that the position 
(x t i p ,2/tip) an d orientation $ of the tip during calcula- 
tions in the laboratory FoR are defined as 

u ih \xt ip (t),yti P {t),t) = u», (17a) 
u (h \x tip (t),y tip (t),t) = u„ (17b) 
*(t) =arg((d x +id y )u {h \x tip (t),y tip (t),t)) (17c) 



and the comoving FoR is chosen so that in it, 
(xtip, 2/tip) = (0,0) and $ = n/2 at all times. Unlike 
other equations of system (16), these are not prescrip- 
tive and is essentially an arbitrary choice, dictated by 
properties of particular systems. We shall refer to the 
pinning conditions (16b, 16c) as "Choice 1", as below we 
shall consider a variation of these, which we call "Choice 
2". 

When h = 0, the system (16) decouples, as its upper 
part including (16a), (16b) and (16c) becomes indepen- 
dent of the lower part (16d). This is the "skew-product" 
decomposition, the upper part describing the dynamics 
in the space of group orbits, so called "quotient system" , 
and the lower part the "symmetry group extension", i.e. 
dynamics along the group, which depends on but does not 
affect the quotient dynamic. The connection between the 
quotient system and the group extension is via the dy- 
namic variables (c, oj); in the following, we refer to these 
three quantities as "quotient data" for brevity. 

The skew-product representation has been useful for 
the analysis of various types of meander of spiral waves 
[29, 30, 37, 38]. Note that the approach used in [30, 38] 
(also see references therein) is based on the assumption 
that the meandering pattern in question is considered 
in the vicinity of a bifurcation from the rigidly rotat- 
ing spiral wave solution, so that the quotient dynamics 
can be reduced to the centre manifold, hence instead of 
equations (16a), (16b) and (16c), these studies consid- 
ered normal forms on the corresponding centre manifolds. 
However, as noted in [39], the Centre Manfold Theorem 
is not applicable for spiral waves, so this approach seems 
to be fundamentally flawed. This technical difficulty of 
course does not in any way affect the validity of system 
(16), which, as we have just demonstrated, is derived by 
elementary means without recourse to any bifurcations. 

In the rest of the paper, we consider system (16) as a 
computational tool, rather than an instrument of theo- 
retical analysis. The disadvantage of the original system 
(1) as a computational tool is that it requires a big com- 
putational grid to simulate dynamics of a spiral in an 
infinite medium, particularly when the tip of the spiral 
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performs excursions to large distances. This is actually 
not necessary, as the dynamics of the spiral is mostly de- 
termined by the events in some finite vicinity of its tip 
[24]. The system (16) takes advantage of this property so 
that the PDE calculations are done always in some fixed 
vicinity of the spiral wave, whereas the movement of the 
tip is described by the ODE part. 



III. NUMERICAL IMPLEMENTATION 

a. Discretization. We use time discretization with 
constant step At and square spatial grid with step A x , 
covering spatial domain (x,y) £ [— L/2, L/2] 2 , so that 

v((i - i )A x , (j - j )A x , kA t ) ~ v^- 
fi8'*|l = l,...,n), i = 0,...JV X) j = 0,...N y , 

N x = N y =L/A x , 

and the grid coordinates of the origin are 

i = (N x + l)/2, j = (N y + l)/2 

(we only use odd values of N x = N y ). We designate 
the k-th time layer of the numerical solution as V fe = 
| % = 1, . . . N X J = 1, . . . N y ) . We discretize the ODE 

dynamic variables on the same time grid, i.e. R(kA t ) ~ 

R k etc. 

b. Operator splitting. We rewrite equation (16a) in 
the form 

= F[v] +H[v;i?,9] + A[v;c,oj] (18) 

where differential operators J 7 , % and A are defined as 
J[v] = DV 2 v + f(v), (19a) 
H[v; R, 6] = h(v, e^ e • Vv, R + e^ e f, t), (19b) 

A[v,c,uj} = (c- V)v + cj|^ 

= (c x - uy) — + (cy +ujx) — . (19c) 

Let LF, T-L and A be discretizations of J 7 , TL and A. Our 
computations proceed as follows: 



c. Kinetics. As specific examples, we consider two 
models, the FitzHugh-Nagumo model [40, 41]: 



yfc+l/3 = yfc + At fr ( yf 



(20a) 



yk+2/3 = yk+1/3 + At? ^ ^ V fc+1/3 , i? fe , , (20b) 

yk+1 = yk+2/3 (yk+2/3 ~k+l > ^fe+l^j ^ ^ 

Mi,2,3 (V fc+1 ) = 0, M 4 (v k+1 ) > 0, (20d) 

6 fc+1 =© fe + A t o; fc+1 , (20e) 

ti k+1 =k k + A t e^ k+1 c k+1 . (20f) 



f : 



u 








V 





a 1 (u — it 3 /3 — v) 
a(u + (3 — jv) 



and Barkley's model [32, 42] 



u 








V 





c 1 u(l — u) (u — (v + b)/a) 
u — v 



(21) 



(22) 



both with D = 



1 


d. Reaction- diffusion step. The computational 
scheme is designed as an extcntion to the standard 
approach to simulation of spiral waves. Specifically, we 
chose Barkley's EZ-SPIRAL [32, 33, 42] as the starting 
point, and extended it to add the other computational 
steps. So the reaction-diffusion step (20a) is as imple- 
mented in EZ-SPIRAL, with central 5-point difference 
approximation of the Laplacian, without any features 
specific to the Barkley model, such as implicit treatment 
of the kinetic terms, and with appropriate modifications 
when FitzHugh-Nagumo model is used. 

e. Perturbations. Wc consider one particular type of 
nonzero perturbation, the electrophoresis, 

h = E9 x u, 

h = E (cos(0)9 x v(r) - sin(e)0 tf v(r)) , (23) 



where E is a diagonal matrix, E 



|E|| < 1. 



E x 
E 2 

For a reaction-diffusion system this perturbation can de- 
scribe movement of the reagents in response to elec- 
tric field with velocities —Ex and — E<z along the x-axis. 
For E = eD, this perturbation can also approximately 
describe the movement of an axially symmetric scroll 
ring. For a cylindrical system of coordinates (r,0,z): 
x = rcos#, y — rsin#, z — z, the diffusion term has 
the form DV 2 u = D(<9 r 2 + \d r + ^d 2 e + d 2 z )u, which for 
dg = and large r is equivalent to an unperturbed diffu- 
sion term with a 2-dimensional Laplacian in (r, z) plane 
plus a small perturbation ^Dd r u. If the filament of the 
scroll is located at large values of r w 1/e and as the dy- 
namics of the scroll is mostly determined by the events 
near its filament, then 1/r can be approximately replaced 
with e. 

Perturbation (23) violates only rotational symmetry of 
the problem, preserving symmetry with respect to trans- 
lations in space and time. Hence h explicitly depends 
only on O. This limitation is not principal and transla- 
tion symmetry breaking perturbations can be considered 
similarly, in which case h would also explicitly depend on 
X, Y and/or t. We discretize the first spatial derivatives 
in the perturbation term using upwind second-order ac- 
curate differences, and use explicit Euler timestepping. 
In the absence of perturbations, h = 0, the perturbation 
step (20b) is of course omitted and V fe + 2 / 3 = V fc+1 / 3 . 



/. Tip definition and pinning conditions. Discretiza- 
tion of the pinning conditions (16b), (16c), using l\ = Z 3 , 
and the right-side first-order discretization of the x- 
derivative, gives 



v 

M 2 ),k 



-(h), k 
J i + l,ja 



It*, 
It*, 



(24a) 
(24b) 
(24c) 



where («o,Jo) are grid coordinates of the origin. This 
works in principle, but gives rather inaccurate and noisy 
approximations for lo, which get worse for finer discretiza- 
tions. This is typical for numerical differentiation. We 
overcome this by enhancing the spatial discretization 
step, by replacing the condition (24c) with 



~(h),k 



(25) 



where the grid point was chosen some way away 

from the centre point, (*i,ji) = (^o , Jo) + (*inc, Jinc 

). This 

means replacing the third pinning condition (16c) with 



v {h) (f inCl t) 



(26) 



where r inc = (A^c, A x j inc ). Empirically, we have 
found that the length of the displacement n nc should 
be of the order of, but not exceeding, one full wavelength 
of the spiral. 

This revised orientation-pinning condition still does 
not define the position uniquely, as illustrated by fig. 2. 
An extra inequality is required to distinguish between 
different solutions satisfying conditions (24a, 24b, 25). We 
use 



corresponding to 



4h),k . 
v h,h <v * 



(27) 
(28) 



Specifically, we chose l\ = I3 = 1 and l 2 = 2. Condi- 
tions (27) and (28) then mean that the third pinning 
condition (25,26) ensures that the front, rather than 
the back, of the excitation wave passes though the grid 
point So equations (16b, 26, 28), with their dis- 

cretizaions (24,27) are our "Choice 2" pinning conditions. 

The Choice 1 and Choice 2 pinning conditions define 
different RMs and different quotient data c(t),u>(t), for 
the same solution u(r, £). However, the two FoRs they 
define have a common origin and differ only by the ori- 
entation angle. So if (c, oj) are quotient data for Choice 
1 pinning conditions, and (c',u/) are quotient data for 
Choice 2 pinning conditions, then we have 



e 7(*-*72)-/ 



d$/dt, 



(29) 



where $ the tip orientation angle in the Choice 2 co- 
moving FoR, so $ — 7r/2 is angle of one FoR against the 
other. 




A 



FIG. 2: (Color online) Non-uniqueness of the revised tip pin- 
ning condition. 



g. Advection. We use an upwind second-order accu- 
rate approximation of the spatial derivatives in A. The 
steps (20c) and (20d) are done in conjunction with each 
other. The discretization of V' £+1 at the tip pinning 
points, resulting from (20c), is used in the three equa- 



tions (20d) to find the three unknowns 



?fc+l pfk+l 



and 



-,fc+i 



so that the pinning conditions (20d) are always 
satisfied exactly (to the processor precision) after every 
step [52]. 

h. Boundary conditions. Since the boundaries in 
the comoving FoR do not represent any physical real- 
ity but are only a necessity of numerical approximation, 
the results can only be considered to be reliable if they do 
not depend on the boundary conditions. So we use both 
Dirichlet and Neumann bondary conditions and compare 
the results. For Dirichlet conditions, we use boundary 
values of the resting state v r , such that f (v r ) = 0. 

i. Tip trajectory reconstruction. The steps (20e) 
and (20f) are simple first-order implementations of the 
corresponding ODEs. The resulting O is used in calcula- 
tions of the % step when the perturbation is on. Other- 
wise, and R are calculated only for the record. 

j. Some details of software implementation. For sta- 
bility purposes, we ensure that the following inequalitcs 
are observed during computations, 



< 



IS I < 
\u\ < 



2A t ; 
A 2 
2 A t ' 
1 



N x A t 



This is an empirical choice motivated by von Neumann 
stability analysis. 

When the absolute values of c x and c y found in (20c) 
and (20d) are beyond these limits then they are restricted 
to the intervals stated above. Also, we eliminated the 
need to restrict the values of c x and c y to their stabil- 
ity limits by moving the spiral wave solution so that the 
tip of the spiral wave is in the center of the box, using 
the standard EZ-SPIRAL's 'mover' function, which per- 
forms translation of the solution by an integer number 
of grid steps, suitably extrapolating the solution where 
necessary near the boundaries. 




For ui, we implemented the restriction that if |w| ex- 
ceeded its maximum stability value, then then we set 
10 = 0. Effectively this means that unless the orientation 
of the spiral wave is already very near the standard ori- 
entation satisfying equation (25) and inequality (27), the 
code computes a solution of the problem 



dv 



— = DV 2 v + f (v) + h (v, Vv, R + e^f, t 

+ (c- V)v, 

w (il) (0,t) = u„ v {h \6,t) = u„, 

dR 
~dt 



(30a) 
(30b) 

(30c) 



instead of (16). That is, it performs reduction by the 
subgroup of translations of the Euclidean group. 

A typical run of the program in the interactive mode 
starts from obtaining a spiral wave solution in the stan- 
dard "ride-off" mode, by solving initial-value problem 
(1). When the spiral wave is initiated so there is one 
tip in the solution, the user switches the program to the 
"ride-on" mode, with calculations according to the above 
scheme. On the switch, the program first of all moves the 
tip of the spiral to the centre of the box via EZ-SPIRAL's 
'mover' function, i.e. by parallel translations of the solu- 
tion, supplementing the missing pieces near boundaries 
by duplicating the existing boundary values. From then 
on, the spiral continues to rotate with its tip fixed at the 
centre of the box, thus solving the problem (30). In this 
regime, only the first two pinning conditions are satisfied, 
and only c x and c y are calculated and used, where as to 
is calculated but replaced with zero, until it falls within 
the stability limit and the fourth inequality-type pinning 
condition is satisfied. From that point, the program pro- 
ceeds in the fully engaged mode, calculating the problem 
(16). 



IV. PRIMARY EXAMPLES: RIGIDLY 
ROTATION AND MEANDER 

First we illustrate how our approach works using two 
examples. One example uses Barklcy model with rigidly 
rotating spiral waves, and the other is FitzHugh-Nagumo 
model with meandering spiral waves. 

Figure 3 illustrates the work of EZRide in the case 
of a rigidly rotating spiral wave. The panels represent 
three consecutive runs, in different regimes: the "direct 
numerical simulations" (DNS) of system (1), then the 
"skew-product" calculation in the comoving FoR, and 
then again the DNS in the laboratory FoR. The skew- 
product calculation in turn consists of two parts. The 
first part is described by (30) where only the two transla- 
tion pinning conditions are engaged, so that the position 
of the tip of the spiral is fixed, but not its orientation, so 
the FoR is co-translating but not co-rotating. The sec- 
ond part is where all four pinning conditions are engaged, 
and the FoR is co-translating and co-rotating. It is seen 




FIG. 3: (Color online) Three consecutive runs of Barkley 
model, a = 0.52, b = 0.05, c = 0.02, L = 20, A x = 1/5, 
A t = 1/2000, r inc = (2,0). The runs t € [0,11] and 
t £ [22,33] are direct simulations. The run t £ [11,22] is 
a quotient system simulation, the pinning points are indi- 
cated by small white crosses. The third pinning condition is 
engaged at t ~ 16.5. 



from fig. 3, that after a transient period, the solution in 
the fully comoving FoR becomes stationary. This corre- 
sponds to the definition of a rigidly rotating spiral wave 
M S ir l relative equilibrium. 

Figure 4 shows a similar set of runs for a different case, 
where the spiral wave is not stationary but is meandering. 
In this case, the solution in the comoving FoR is not 
stationary, but periodic in time. This corresponds to 
the definition of a meandering spiral wave as a relative 
periodic orbit. 

Figure 5(a,b) show selected pieces of tip trajectories 
obtained as a result of the runs shown in fig. 3 and fig. 4. 
The discretization steps there are deliberately chosen 
crude, to allow very fast running simulations, and also 
to illustrate the difference introduced by the change of 
method of computation. The tip trajectories obtained by 
reconstruction from the quotient data are qualitatively 
similar to the tip trajectories obtained in DNS. However, 
the quantitative difference is also quite evident. In the 
case of rigid rotation, the reconstructed trajectory radius 
is noticeably bigger than that from DNS, and the centres 
of the meandering patterns in different runs arc offset 
against each other. As panels (c,d) in the same figure 
show, these discrepancies decrease when the discretiza- 
tion steps are refined. 

Figure 6 shows the tip and quotient system trajecto- 
ries, obtained in laboratory and comoving FoR calcula- 
tions, for a meandering spiral. This is drawn for the 
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FIG. 4: (Color online) Three consecutive runs of FHN model, 
a = 0.2, p = 0.7, 7 = 0.5, L = 30, A r = 1/3, A t = 1/720, 
r inc = (20/3,0). The runs t £ [0,22] and t € [44,66] are 
direct simulations. The run t £ [22, 44] is a quotient system 
simulation, the pinning points are indicated by small white 
crosses. The third pinning condition is engaged at t m 27.5. 

finer discretization steps, as in fig. 5(d). For comparison, 
quotient data for both the laboratory and comoving FoR 
calculations were recalculated for the Choice 1 pinning 
conditions using (29). There is good agreement between 
the two methods of calculations, within the expected ac- 
curacy. More detailed analysis of the numerical accuracy 
of our method is given in the next session. 



V. NUMERICAL CONVERGENCE 

Figure 7 illustrates the convergence of the results of 
calculations of rigidly rotating spiral, using EZRide with 
Neumann and Dirichlet boundary conditions, and DNS 
using Neumann boundary conditions. In these calcula- 
tions, the box size is fixed at L = 60 and the timestep 
is changed with the spacestep so that A t = A^/40. For 
w(A^) dependence, we also show the angular velocity 
measured in direct numerical simulations. We do not 
show |c(A^)| found in DNS, since obtaining it involves 
numerical differentiation which gives accuracy insuffi- 
cient for the convergence study. 

Our discretizations are second order accurate in A x 
and first order accurate in At both in DNS and in the 
riding mode, which corresponds to linear dependence of 
any results on A x for A x —> 0. We see in fig. 7 that this is 
indeed the case. Linear extrapolation of the w(A|) gives 
the values of w(0) for laboratory and comoving calcula- 
tions coinciding to within 10 -3 . 




i i i i I " L i i --f i 

-4 -2 2 4 -2 -1 1 

(C) X (d) X 

FIG. 5: (Color online) (a,b) Reconstructed tip trajectories 
from (a) simulation shown in fig. 3 and (b) simulation shown 
in fig. 4. The pieces labelled 1 are trajectories obtained in di- 
rect simulations in the laboratory FoR. The pieces labelled 2 
are trajectories obtained via co-translating simulations, with 
first two pinning conditions engaged. The pieces labelled 3 
correspond to co-moving (co-translating and co-rotating) sim- 
ulations with all three pinning conditions engaged. The final 
pieces labelled 4 correspond to direct simulations in a non- 
moving FoR, which has been displaced with respect to the 
laboratory FoR during the quotient system simulations, (c) 
Same as (a), with A x = 1/10, At = 1/4000. (d) Same as (b), 
with A x = 1/10, At = 1/4000. 



One of the advantages of EZRide is the fact that the 
simulations can be done in a smaller box compared to 
DNS. So, the last test is convergence in box size. We 
have calculated the rigidly rotating spiral by EZRide at 
fixed A^ = 1/15, A t = 1/9000 and L varying through 
[15,60] and found that both |c| and |a>| vary by less than 
10-3. 



VI. APPLICATION I: THE 1:1 RESONANCE IN 
MEANDERING SPIRAL WAVES 

One of the cases where the DNS would meet with diffi- 
culties, is the study of the the meandering of spiral waves 
for parameters near the "1:1 resonance" between the Eu- 
clidean and the Hopf frequencies. This case is marginal 
between meandering patterns with inward petals and 
outward petals. Near the resonance, the spatial extent 
of the meandering trajectory becomes large, and for the 
case of exact resonance, infinite, and the spiral appears 
to be spontaneously drifting [25, 43]. Hence, following 
the dynamics of the spiral wave in the comoving FoR 
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FIG. 6: (Color online) Meander in the FHN model, calculated 
in the laboratory frame of reference (DNS), and from quotient 
system (EZRide). In (a), the meandering pattern is shown, 
which for the EZRide curve is obtained by numerical integra- 
tion of quotient data using (15). In (b-d), the projections of 
the quotient data are shown, which for the DNS curves are ob- 
tained by numerical differentiation of the tip trajectory, using 
(15). 
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FIG. 7: (Color online) Convergence of the rigidly rotating 
spiral wave solution in the Barkley's model. 



presents an advantage. 

We illustrate this using the FHN model. We fix the 
discretization parameters at A x — 1/8, A t = 1/2560 and 
L = 20. The choice of model parameter is influenced 
by Winfree's "Flower Garden" [44], which gives a rough 
estimate for the location of the 1:1 resonance line in the 
(a, /3) plane at 7 = 0.5. Using this information, we have 
selected two values a = 0.2 and a = 0.25, and scanned 
values of P across the resonance value, which we deter- 
mined as fa « 0.93535 for a = 0.2, and (3 ~ 0.81362 for 
a = 0.25 at our discretization parameters. 

The results are presented on figures 8-11. The shape of 
trajectories is well known from the theory, and is outward 
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FIG. 8: (Color online) The reconstructed tip trajectories in 
FitzHugh-Nagumo system with a — 0.2, 7 = 0.5 and varying 
P- 



petals for P < p and inward petals for P > p , degener- 
ating into spontaneous straightforward drift at P = fto- 
The trajectory at P = Pq in fig. 8 is shown twice: once 
for the whole duration as it was calculated, fig. 8(c), and 
then a close-up of small part of it, fig. 8(d). Calculation 
of this particular trajectory using DNS would require, by 
our estimate, about five weeks, as opposed to 2.5 hours 
used by EZRide. 

The change of the quotient dynamics with changing P 
is illustrated in fig. 9. As opposed to the tip trajectories, 
there is no evident qualitative changes in the shape of 
the limit cycle across P = Pq. Note the very elongated 
shape of the limit cycles in all three projections. We do 
not know whether this has some theoretical explanation 
or is merely incidental. 

The parametric line a = 0.25 exhibits similar be- 
haviour, as shown in fig. 10 and 11. This is closer to 
the Hopf bifurcation line in the quotient system, called 
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FIG. 9: (Color online) Various projection of the limit cycles in the quotient system corresponding to the trajectories shown in 
fig. 8. 



dM line in [44]. Correspondingly, the size of the limit 
cycles in the quotient system is smaller and they become 
more oval-shaped Note that the scale of c y -axis is dis- 
proportionately stretched in fig. 11, i.e. the Hopf central 
manifold appears to be nearly orthogonal to that axis. 
Again, there is no qualitative change in the quotient sys- 
tem dynamics when crossing the 1:1 resonance. 

VII. APPLICATION II: LARGE CORE SPIRALS 

Another example where the spatial extent of the spiral 
wave dynamics is large is the vicinity of Winfrec's "ro- 
tors boundary" dR in the parametric space [44]. In the 
vicinity of this boundary, the period of rotation and the 
radius of the core of the spiral wave grow infinitely. 

There are at least two different asymptotic theories, 
based on different choice of small parameters, which aim 
to describe the vicinity of dR. Hakim and Karma [45, 
46] have developed a "free-boundary" asymptotic theory 
applicable to FitzHugh-Nagumo type models in the limit 
c — > or a — > in terms of our chosen kinetics, where 
angular velocity to typically decreases as 

3 

\u>\ oc \p-p+\?, p-^p*, (31) 

where p is a parameter of the model such that p = p* 
corresponds to the dR boundary. 

Elkin et al. [47] obtained an alternative asymptotic 
based on assumptions which were not restricted to kinet- 
ics of any particular kind, but which were not directly 
validated. Their prediction was 

\w\ oc \p-p*\, p^p*. (32) 

Further analysis has suggested that these two alter- 
natives are not actually antagonistic and may be even 
observed in the same system in different parametric 
regions [48]. Reliably distinguishing between the two 
asymptotics is challenging for DNS as it requires a rather 
close approach to the critical point p = p t , which is not 
known a priori, implying large tip trajectory radii and 
correspondingly significant computational resources. 



In here we present an example of studying this depen- 
dence using calculations in the comoving FoR, which is 
free from the above complication, as it can be performed 
within the box of fixed size for all p. 

For this study, we use Barkley's model with varying 
parameter p chosen to be a, varying from a = 0.48 down- 
wards with step 0.001 until 0.43, with other parameters 
fixed at b = 0.05 and c — 0.02. The discretization pa- 
rameters are L = 30, A x = 1/8, A f = 1/2560 and 
Fine = (0,7/4). 

Selected stationary solution obtained in this way are 
illustrated in fig. 12, and the graphs of cu(a) and c y (a) 
are shown in fig. 13. We compare the features of the 
observed solutions with those that are given by the two 
asymptotic theories [48], and observe that 

1. There is a critical value of the parameter a* « 
0.456, at which the behaviour of the solution 
changes qualitatively. At a = a*, we observe a 
nearly straight broken excitation wave. 

2. For a > a», the solutions arc spiral waves, that 
is, broken excitation wavelets, which become less 
and less convex as a — > a*, and have macroscopic 
angular velocity which however diminishes in the 
same limit. 

3. For a < a* the solutions are retracting nearly 
straight but slightly concave wavelets, with very 
small angular velocity. 

4. For a = a*, the direction of movement of the tip 
seems approximately orthogonal to the overall ori- 
entation of the wave itself. 

5. For a < a*, the vertical component of vector c de- 
pends on a in a way which is consistent with the 
asymptotic \c y \ oc |a — a*! 1 / 2 , see fig. 13(e,f). Since 
the overall orientation of the wavelets, as seen in 
fig. 12(a-c), is nearly vertical we can take c y as a 
crude estimate of the "global tip growth rate" as 
defined in [48]. 

6. For a > a», the angular velocity of solutions de- 
pends on a in a way which is consistent with the 
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FIG. 10: (Color online) Same as fig. 8, for a = 0.25. 



asymptotic \u>\ oc a — a*|, see fig. 13(b) but not 
|w| oc |a — a*| 3 / 2 , see fig. 13(c). 

All these observations are in agreement with the the- 
ory in [48] and can be used to empirically distinguish be- 
tween the Elkin et al. asymptotics (corresponding to the 
"I/V" parametric boundary in [48]) and Hakim-Karma 
asymptotics (respectively, "J/C" boundary in [48]). 

Feature 1 is inconclusive: existence of a critical so- 
lution, called "critical finger" by Hakim and Karma, is 
common for both J/C and I/V boundaries, but the shape 
of this solution is different. It is asymptotically linear for 
I/V boundary, and asymptotically logarithmic for J/C 
boundary. Looking at fig. 12(d) and considering the ef- 
fect of the boundary conditions, it is not clear which case 
is nearer to the observed reality. 

Feature 2 is common for I/V and J/C boundaries. The 
phenomenological difference is that spirals close to I/V 
boundary can be "growing" or "shrinking" , while spi- 
rals close to J/C boundary can only be "growing" . The 



movement of the tip in fig. 12(d-f) seems approximately 
orthogonal to the orientation of the wavelet near the tip, 
which is consistent with both cases. 

Feature 3 tips the balance in favour the I/V bound- 
ary since the broken wavelets are concave. According to 
[48] , the translating waves near an I /V boundary should 
be concave, and those near an J/C boundary should be 
convex. 

Feature 4 is common for I/V and J/C, as in both cases 
the critical fingers should have zero "global growth rate" . 

Feature 5 is common for I/V and J/C boundaries. 

Feature 6 is, in our opinion, a convincing evidence in 
favour of an I /V boundary, since according to [48] , near 
I/V boundary the dependence lo(8) is linear, whereas 
near J/C boundary it is |u;(5)| oc |<5| 3 / 2 . 

An unequivocal interpretation of all theoretical pre- 
dictions in the view of our present numerical results 
would require further investigation, as the asymptotics 
of [47, 48] operate with a "crest line" of an excitation 
wave. There is no obvious operational definition of this 
line which would be valid up to the tip, and some of 
the predictions concern the mutual orientation of this 
line and the tip velocity. However the predictions that 
can be unambiguously interpreted, seem to indicate that 
for the model considered here, we have the case of I/V 
boundary, i.e. Elkin et al. asymptotics, rather than J/C 
boundary corresponding to Hakim-Karma asymptotics. 

The last observation here is that of the small angular 
velocity uj calculated for the "retracting waves" at a < a* , 
seen on fig. 12(a-c). As we already noted, the small- 
ness of these w values is consistent with the theoretical 
prediction of translating but not rotating waves. How- 
ever when these values are magnified, we observe that 
they demonstrate a peculiar power law |w(a)| oc |a — a„\ p 
where p ~ 1/4.3, see fig. 13(d). A theoretical explanation 
of this requires further study; it is clear, however, that co 
in this area is strongly affected by the boundaries, as the 
curves for L = 30 and L = 35 differ quite significantly. 



VIII. APPLICATION III: ELECTROPHORESIS 
OF MEANDERING SPIRAL 

Finally, we illustrate calculation of the movement of 
spiral waves in a perturbed reaction-diffusion system. We 
consider FitzHugh-Nagumo kinetics at the same param- 
eters as in fig. 4, and add to it the "electrophoresis" per- 
turbation (23) in the right-hand side, with E = eD. 

Results of the simulations are presented in fig. 14. The 
unpertubed spiral waves for these parameters are me- 
andering, so with the perturbation present, we observe 
meandering with drift. The drift proceeds with a con- 
stant average velocity, which is consistent with the fact 
that the perturbation violates only the rotational but not 
the translational symmetry of the problem. The average 
drift is to the left, which corresponds to a collapsing scroll 
wave. So at these parameter values, the scroll waves have 
positive tension, inasmuch as this concept can be applied 



12 





a = 0.450 



1. 



a = 0.460 



FIG. 12: (Color online) Snapshots of relative equilibria in 
Barkley model obtained at different values of parameter a. 
The arrows indicate the direction of the vector c. 



to meandering scrolls. 

In the calculations in the laboratory FoR, the time 
during which the drift can be observed is limited, as when 
the spiral reaches the left boundary, it terminates. In the 
comoving FoR, this drift can be observed indefinitely. 
Comparing the traces in fig. 14 we see that although, as 
we know from figures 5 and 6, the discretization is too 
crude to give quantitative agreement between laboratory 
and comoving calculations in detail, the drift velocities 
obtained in these two ways are very similar. 

We illustrate the relative advantages of the two meth- 
ods of calculation by comparing the computation costs. 
The laboratory FoR simulation, for L = 30 and t G 
[0, 300] has taken 325 sec (the spiral has annihilated at 
the left at t ~ 237). The time taken by the comoving FoR 
simulation for the same boxsize L and the same t interval 
is 462 sec, i.e. is naturally somewhat longer due to the 
extra effort required for the advection term calculations. 
However, the comoving FoR calculation proceeded un- 
abated where the laboratory FoR calculation failed due 
to annihilation with the border. To continue the labora- 
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FIG. 13: (Color online) Dependencies u>(a) and c y (a) of the 
relative equilibria, for different L as indicated. On panel (a), 
the symbols correspond to the selected values of a used in 
fig. 12. 



tory FoR calculation to the same extent we would have to 
increase the box size L with a corresponding increase in 
computation cost. Moreover, virtually the same result, 
as far as drift velocity is concerned, can be obtained by 
comoving FoR calculation with L = 20, and it takes only 
202 sec. Of course the drift in the laboratory FoR with 
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FIG. 14: (Color online) Trajectories of tips of drifting meandering spirals calculated in the laboratory FoR (for L = 30) and 
in the comoving FoR (for L = 30 and L = 20). The thin black dotted lines designate the boundaries of the calculation box 
in the laboratory FoR where the initial position of the tip is in the centre. The parameters are the same as in fig. 4 and the 
perturbation is h = eDc^u, where e = 0.1. 



L = 20 would terminate even earlier. 



IX. DISCUSSION 

We have described a numerical method of solving a 
reaction-diffusion system of equations describing a spiral 
wave, in a frame of reference which is moving with the 
tip of that wave. 

We have shown the method can provide accurate solu- 
tions, and that there are applications where the computa- 
tional cost of our method can be considerably lower than 
that of the conventional approach, or the conventional 
approach is just inapplicable. As always, the compu- 
tational advantages are particularly essential in case of 
parametric studies, for which the method is well suited. 

Although the applications were chosen just to provide 
some meaningful examples of use of the method, the re- 
sults obtained there can be of scientific value themselves. 

So, we have investigated the vicinity of the "1:1 res- 
onance" manifold in the parametric space, which corre- 
sponds to spontaneous drift of spirals, and which sepa- 
rates meandering patterns with outward petals and in- 
ward petals. Henry [49] has proposed a theory which 
implies that this manifold coincides or is an analytical 
continuation of the manifold where the filament tension 
of scroll waves vanishes. There are reports in literature 
confirming that change of sign of filament tension is as- 
sociated with change from outward to inward petals in 
meandering patterns, but also examples where there arc 
no such correlation, e.g. [50] and references therein. Our 
simulations indicate that as far as orbit manifold dynam- 



ics of the spiral is concerned, the 1:1 resonance is not 
characterized by any special features. Hence any special 
features of this resonance ought to be due to the Eu- 
clidean extension of the orbifold dynamics. Since scroll 
filament tension can also be defined via properties of the 
spiral wave solutions within the comoving FoR, any ge- 
netic and generic relationship between the two manifolds 
seems unlikely (but, of course, cannot exclude the possi- 
bility of such relationship in some special cases). 

We have also investigated the vicinity of the "dR v 
manifold in the parametric space, which has provided 
a strong evidence towards one of the two theoretical pos- 
sible asymptotics, namely Elkin et al. [47] asymptotics 
as opposed to Hakim-Karma [45] asymptotics. It should 
be noted here that while Hakim-Karma asymptotic the- 
ory was based on assumptions which have been well es- 
tablished, the Elkin et al. asymptotic theory was using 
assumptions, validity of which could not be asserted at 
that moment. Here we have presented firm evidence that 
Elkin et al. asymptotis is not a mere theoretical possibil- 
ity but is indeed observed in reality (see also [36] and 
a discussion below). A direct confirmation would be via 
calculation of the "response functions" , i.e. critical eigen- 
functions of the adjoint linearized operator of the critical 
finger solution. This would require obtaining first a good 
quality critical finger solution, so the method described 
here can be a significant step towards this goal, too. 

Finally, we have demonstrated that calculations in 
the comoving FoR can be efficiently used to study 
perturbation-caused drift of spirals, including meander- 
ing spirals. Although the asymptotic theory of drift of 
meandering spirals is yet to be developed (see, however, 
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a preliminary draft of such theory in [51]), we can expect, 
for instance, that scroll waves in the FitzHugh-Nagumo 
model with the parameters as in fig. 4, 14 will have "pos- 
itive tension", i.e. tend to collapse, rather than develop 
a scroll wave turbulence. The advantage of calculating 
drift in the comoving FoR, apart from computation cost, 
is absence of "pinning" effects of spatial discretization, 
both in terms of discrete space steps and discrete spatial 
directions, on the drift. 

Our approach can be compared to the approach pro- 
posed by Beyn and Thummler (BT) [35] . BT use a sim- 
ilar mathematical idea of decomposing the evolution of 
the nonlinear wave into the motion of the wave and evo- 
lution of its shape, which in the functional space appears 
as decomposition into motion along and across the Eu- 
clidean symmetry group orbits. But there are also dif- 
ferences. There are technical details of implementations 
which arc probably of lesser importance, such as choice 
of polar vs Cartesian grid, central vs upwind discretiza- 
tion of spatial derivatives and explicit vs semi-implicit 
discretization in time. More significant differences are 
in the "phase conditions" they use, which play the same 
role as, but are qualitatively different in nature from, our 
"pinning conditions" . One aspect is that the phase con- 
ditions involve integral functionals. We show here that 
this is not necessary, and local conditions like (24) are 
simpler. The other aspect is the one we discussed in 
section II: the BT phase conditions appear to be well 
suited for calculation of relative equilibria (rigidly rotat- 
ing spirals) but not necessarily for relative periodic solu- 
tions (meandering spirals). Further, the phase conditions 
proposed by BT were not intended for use with symme- 
try breaking perturbations that produce drift of spirals. 
And indeed, BT comment in their paper that "it seems 
quite a challenging task to freeze drifting spirals or recog- 



nize meandering spirals as periodic orbits." As we have 
demonstrated, our approach works both for meandering 
spirals and for drifting spirals. 

After completing this study we became aware of a work 
by Hermann and Gottwald (HG) [36] who also investi- 
gated the large core limit, using a further development 
of the BT method. HG have paid a great deal of attention 
to refining the boundary conditions so as to minimize the 
effect of boundaries onto the quotient dynamics. This has 
allowed them, in particular, to verify the linear scaling 
law (32) for seven decades of variation of |w|, compared 
to mere one decade as shown in fig. 13. Notice that as 
shown in the same figure, our progress towards smaller 
values of \ui\ is limited precisely by the influence of bound- 
aries. HG also have explicitly addressed the issue of the 
numerical stability of the computations, which we treat 
in this study purely empirically. 

We believe that combining the advantageous features 
of the approach developed by BT and HG, and the one 
proposed here, is an interesting topic for future work, 
which may yield further results about spiral wave dy- 
namics, that are not possible, or very difficult, to obtain 
by direct numerical simulations. 
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